#pragma once
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp>
#include <glm/gtc/type_ptr.hpp>
#include <glm/gtx/rotate_vector.hpp>
#define _USE_MATH_DEFINES
#include <math.h>

class CubicKernel
{
	protected:
		static double m_radius;
		static double m_k;
		static double m_l;
		static double m_W_zero;
	public:
		static double getRadius() { return m_radius; }
		static void setRadius(double val)
		{
			m_radius = val;
			static const double pi = static_cast<double>(M_PI);

			const double h3 = m_radius * m_radius * m_radius;
			m_k = static_cast<double>(8.0) / (pi * h3);
			m_l = static_cast<double>(48.0) / (pi * h3);
			m_W_zero = W(glm::dvec3(0.0, 0.0, 0.0));
		}

	public:
		//static unsigned int counter;
		static double W(const glm::dvec3& r)
		{
			//counter++;
			double res = 0.0;
			const double rl = glm::length(r) ;
			const double q = rl / m_radius;
#ifndef NO_DISTANCE_TEST
			if (q <= 1.0)
#endif
			{
				if (q <= 0.5)
				{
					const double q2 = q * q;
					const double q3 = q2 * q;
					res = m_k * (static_cast<double>(6.0) * q3 - static_cast<double>(6.0) * q2 + static_cast<double>(1.0));
				}
				else
				{
					res = m_k * (static_cast<double>(2.0) * pow(static_cast<double>(1.0) - q, 3));
				}
			}
			return res;
		}

		static glm::dvec3 gradW(const glm::dvec3& r)
		{
			glm::dvec3 res;
			const double rl = glm::length(r);
			const double q = rl / m_radius;
#ifndef NO_DISTANCE_TEST
			if (q <= 1.0)
#endif
			{
				if (rl > 1.0e-6)
				{
					const glm::dvec3 gradq =  ((double)1.0 / (rl * m_radius))* r;
					if (q <= 0.5)
					{
						res = (double)(m_l * q * (3.0 * q - 2.0)) * gradq;
					}
					else
					{
						const double factor = static_cast<double>(1.0) - q;
						res = (double)(m_l * (-factor * factor)) * gradq;
					}
				}
			}
#ifndef NO_DISTANCE_TEST
			else
				glm::vec3(0,0,0);
#endif

			return res;
		}

		static double W_zero()
		{
			return m_W_zero;
		}

};

